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Introduction 

This work is motivated by NASA plans to develop and maturate technology for autonomous 
cryogenic propellants management (CPM) on the ground and in space. The latter can be defined 
as the capability to accomplish launch vehicle propellant load and drain without human interaction. 
Importantly, the system and subsystems should be able to change their behavior in response to un- 
anticipated events. The solution of this problem requires a deployment of an intelligent control 
system that can optimize loading protocol, predict various deviations from the nominal regime, 
detect, localize, and mitigate faults online. 

The difficulties in solving this problem stem from the complexity of the transient analysis of the 
two-phase flow (TPF) under strongly non-equilibrium conditions. The application of such analysis to 
cryogenic propellant management poses additional complications due to uncertainties of flow 
regime and heat transfer characteristics of the cryogenic fluids. These uncertainties are even larger 
under microgravity condition imposing severe design concern for the development of the 
autonomous CPM in space. 

To solve this problem and to enable online autonomous control of cryogenic loading we have 
developed a concept of physics based automated loading operation in space and on the ground 
that combines hierarchy of models with best industrial practice and deep analytical insight into the 
nature of flow patterns, frictional losses and heat transfer correlations. The model hierarchy is 
ranging from the fastest online fault detection, isolation, and mitigation (FDI&M) model to the full 3D 
supercomputer model on the ground and incorporates (worth of 50 years R&D) best practice of 
nuclear reactor and refrigerator industrial modeling frameworks for the automated control of the 
two-phase flow and heat transfer. 

Flere we report the progress in development of one of the key component of the model 
hierarchy, namely one-dimensional model of fully separated non-equilibrium two-phase cryogenic 
flow that allows for real time control of transient response of the system during chilldown and 
loading. 

In developing this model, we took into account additional limitations imposed on the program by 
constraint resources and time. We recognized that the development of multiphase cryogenic flow 
models capable of predicting flow regimes and heat transfer at low gravity is one of the top priority 
recommendations made by NRC Committee [1] for successful NASA Space Program. In these 
recommendations, it was also pointed out [1] that such models should incorporate best practice of 
autonomous thermal-hydraulic management developed by the Department of Energy. 

Indeed, cryogenic chilldown and loading shares many common transients with other industrial 
processes [2] such as loss of coolant and other anticipated transients without scram 1 . Many 
transient characteristics of system failure during cryogenic loading are also common with industrial 
processes including, e.g. loss of coolant, loss of feeding fluid, loss of offsite power, station blackout 
etc. We note in addition that the latest in the series of nuclear reactor codes are highly generic and 
can be used for simulation and control of a wide variety of hydraulic and thermal transients in both 
nuclear and nonnuclear systems involving mixtures of vapor, liquid, non-condensable gases, and 
nonvolatile solute. 

In what follows, we briefly formulate a general transient two-phase flow problem and list the key 
difficulties in its solution. Next, we briefly discuss the hierarchy of models that can be used to 
accomplish online control of this flow in various approximations. 


the failure of the emergency shutdown system 
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Two-phase flow problem 


At present, the predictive capability and physical understanding of the two-phase flow in 
industrial systems must rely heavily on theoretical and computational models [3], Major theoretical 
efforts here are required to find functional form of correlations that cover the diversity of multiphase 
flows patterns, friction and heat transfer characteristics. 

For autonomous control the model has to be efficient enough to detect and localize faults in real 
time. We therefore restrict our discussion to one-dimensional models that satisfy this requirement. 
We note though that the models discussed below can be extended to 3D version if required (see 
e.g. [2], [4], [5]). 

We note that the hierarchy of one-dimensional models of two-phase flow is diverse. More 
importantly, an efficient online fault detection, isolation, and mitigation must involve optimal 
utilization of a few members of this hierarchy for every particular fault. For example, the most 
accurate model (separated non-equilibrium two-phase flow) at the top of this hierarchy should be 
responsible for real time FDI and providing recommendations for the most likely cause of the fault 
and possible mitigation strategies, therefore restricting search space for the fault and control 
parameters. The fastest model at the bottom of this hierarchy (incompressible isothermal flow 
network) should be responsible for providing guess for the initial conditions and fast estimations of 
pressure losses and flow rates. The intermediate level model (e.g. moving front model for 
homogeneous flow) performs evaluation of relatively slow thermal fault transients and search for 
the optimal mitigation strategies. 


wall T w 



Figure 1. Control volume of the two-phase flow. cc g [ void fractions for gas (g) and liquid (l). T g \ phasic temperatures. phasic 
velocities. wetted perimeter for each phase. I; perimeter of the interface. A is the cross-sectional area. T w is the wall temperature. 

The relation between the models in the hierarchy their capabilities and the conceptual 
architecture of the intelligent online control module should become more transparent in the 
following sections after we describe generic two-phase flow models and basic solution algorithms. 
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Generic two-phase flow model 


The governing equations for effective-field formulation of multi-fluid/multi-phase flows were 
developed in late 60s [4], [6]-[9][6], and are widely used in nuclear reactor safety and thermal- 
hydraulics analysis, as well as for other multiphase-flow applications in oil and aerospace 
industries, geoscience, advanced weaponry, etc.. 

Consider two-phase flow equations averaged over control volume [6], [10], [1 1] shown in the 


Figure 1 . Each phase of the fluid is characterized by its own void fraction a k , density p k , 
temperature T k (energy e k ), and velocity u k , where index k takes values g for the gas/vapor phase 
and l for liquid phase. For non-homogeneous (u a ^ U;) and non-equilibrium (T a =£ T;) flow 
conditions a closed system of equations can be obtained assuming equal local pressure values for 
the two phases (p a = p; = p) that the source terms on the right-hand sides of the balance equations 
are exclusive algebraic functions of state and flow parameters [6]-[8], [10], [11]. The corresponding 
model often referred as “Wallis model” [7] has the following form 


where the first three equations represent mass, momentum and energy conservation laws for 
gas/vapor and the last three equations represent the same laws for the liquid. 

This set of equations can be closed by adding volume conservation equation 


In addition, the constitutive equations for the interfacial mass flux r CT , friction and heat transfer at 
the wall and at the interface have to be given for various flow regimes. The functional form and 
parameters of these constitutive relations are the major unknowns of the two-phase flow problem in 
any specific environment as will be discussed in more details later. 

The main applications of this model are real-time fault detection and isolation, parameter 
estimation, off-line analysis of the loading regimes, analysis and prioritization of the faults. 


To obtain the next model in the hierarchy we add together momenta equations for each phase 
and assuming equal phase velocities (u a = u t = u). Using volume conservation equation and 
introducing mixture density 



( 1 ) 


a g + a, = 1 


and phasic state relationships 



Two-phase flow model hierarchy in one dimension 


Homogenous non-equilibrium model 
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Pm= a S P S +a ,Pl 


we obtain the following mixture momentum equation 

( A Pm U ),, + ( A py), X + A P, X = ~ A P m § Z , X ~{ T X) 2 t ( 2 ) 


In deriving this equation, it was assumed that all the interfacial terms cancel each other giving 
rise to so-called jump condition 

~ T g J, + Ar g \ + - AT g u u = 0 • 

In addition, it is usually assumed without derivation that the interface momentum transfer due to 
friction and due to mass transfer independently sum to zero, that is, 

T gl =T u =T l > U ig =U a = Ui . 

The remaining equations can be left unchanged leading to the so-called five equations non- 
equilibrium homogeneous model (see e.g. [8]). The eigenvalue analysis [8] shows that this system 
has three distinct eigenvalues 

/lj 2 = u±a, Ay =u (3) 


where the first two eigenvalues correspond to pressure waves with sound velocity 




a iP m 
a] Pi 


and the last eigenvalue corresponds to void fraction wave. Two more eigenvalues are equal T 3 and 
correspond to the entropy waves. 


Homogeneous equilibrium model 

The homogeneous equilibrium model can be obtained from the homogenous non-equilibrium model 
by assuming that two phases have the same temperature ( T g = 7) = T ), velocities (u a = = u) 

and pressure (p a = p t = p) and by adding pairwise the equations of conservation for the mass and 
energy. The resulting set of equations has the form 

A PmS+( A Pm U )s=°> 

A (py,+( A py) x =- A p, x -yL) 2 t-p>" Agsin0 ’ (4) 

A {Pm E ), t +{ A Pm UH ), x = 


where E = e + u*/2 is total energy, /?=£■ + p/p = h + u 2 /2 is the total enthalpy, & is the angle of pipe 
axis with horizontal line, x is the length along the pipe, and the mixture properties are defined as 
follows 

P = a gP g + a iPn e = X 8 e g +Xi e i’ h = z g h g +zA (5) 

Here z a (v > ' s the corresponding mass fraction 

a g(i)Pgd ) 

/£#(/) ’ 
a g P g +ot,p, 

and are specific energy and enthalpy of gas (liquid). 

In one-component two-phase media undergoing phase transitions (evaporation or condensation), 
the thermal equilibrium assumption results [8] in a strict coupling of pressure and temperature as 
given by the saturation condition p = p saC (T ) as long as 0 < j < 1. The structure of three 
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eigenvalues remains the same as in eq. (3). The primarily application of this model is a fast and 
reasonably accurate solution of the online optimization problem for thermal transients in non- 
isothermal compressible TPF networks. 


Steady flow network 


The last member of the considered hierarchy is the isothermal flow network model. It is 
assumed that the flow in the network has single phase and is incompressible and isothermal. Then 
the last equation in (4) is reduced to the hollowing statement that the fluid temperature in the 
network remains constant 

T = const, 

while the fist equation in (4) is reduced to the conservation of the total mass flow rate and 
volumetric flow rate at every cross-section and junction 

Ap m u = const , p m = const, and Au = const . 


Finally, under assumption of steady flow the momentum equation in (4) integrated over section 
of the pipe with length L and diameter D becomes 

r fL 


Ap = - 


' pu 2 


D 


-pgAh 


This case is important for a few reasons. Firstly, it allows for the fast non-linear optimization during 
the search of the fault and control parameters. Secondly, this is the simplest possible 
approximation of the problem that yields useful results. In addition, this model is currently used by 
the KATE for the online FDI&M at KSC cryogenic testbed (CTB). 


Intelligent autonomous loading control 


An analysis of the current state of art modeling capabilities of the two-phase flow suggests that 
an efficient online autonomous control of the TPF can be accomplished in an optimal way using 
hierarchical approach to the control architecture. Within this architecture the most accurate of the 
available models (separated non-equilibrium TPF model) follows in real time the loading operation, 
detects the deviations from the nominal regime, generate a table of suspected faults and 
corresponding sensitivity matrixes, suggests the most probable mitigation strategies and the 
efficient method of using module capabilities for each specific fault. The most important use of this 
model is the detection of the critical faults at the earliest possible time. 

The primarily use of the low level models at the bottom of the hierarchy (steady compressible 
and incompressible flow networks) is to provide estimations for the flow rates and pressure 
distribution along the system conditions in two limiting cases (compressible and incompressible 
flow). The corresponding estimations and mitigation strategies will be used as initial conditions by 
more accurate models. 

The key application of the middle level models (e.g. homogeneous moving front model) is to 
provide the fastest possible solution of the online of optimization problem for the fault detection, 
evaluation, and mitigation within specified parameter’ and control’ subspaces. Essentially, for the 
cryogenic loading the model has to provide analysis of the temperature transients both for the wall 
and for the fluid and to optimize mitigation. 

A very important application of the control module is off-line analysis and prioritization of each 
deviation of the loading operation from the nominal regime. The module will collect and analyze 
severity of the detected faults, continue to learn model parameters (physics of the two-phase flow) 
in autonomous regime. The module will also determine if further analysis of the loading regimes 
and/or system/components performance is required using ground facilities including high-fidelity 
supercomputer code and/or experimental testbeds. 
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In this report, we focus on the development of ine of the key components in this hierarchy -- 
fully separated non-equilibrium two-phase flow - which will be used for the real-time fault detection, 
isolation and parameter estimation. In the next section, we review briefly the current state of art in 
the methods of numerical solution of this model. Then we will describe briefly one of the solution 
methods, which is currently used for modeling cryogenic loading at KC CTB. 


Solution methods (in a few words) 


In our brief discussion below, we follow closely recent review [4], The reader should consult this 
comprehensive text for further details. All current work-horse reactor thermal-hydraulics codes 
(RELAP5 [2], TRAC [12], TRACE [5], CATHARE [13] and RETRAN [14]) originate from Liles and 
Reed [15] extension of Harlow and Asden [16], [17] all-speed “Implicit Continuous-Fluid Eulerian 
(ICE)” algorithm. The whole family of these algorithms is often referred to as “semi-implicit”. 
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Figure 2 Hierarchy of the ID TPF models and relation of the algorithms (see [4] for further details). 

A couple of extensions of the “semi-implicit” algorithm were introduced, based on the fractional 
step method, to enhance stability of the method and to eliminate material CFL restrictions. These 
are the SETS algorithm (implemented in TRAC and TRACE, [5], [12]), and Nearly-lmplicit algorithm 
(implemented in RELAP5-3D, [2]). 

We refer to ICE-based algorithms as “weakly-compressible”, even though the compressibility 
effects are fully accounted for in numerical discretization. Almost all reactor thermal-hydraulics 
codes belong to this class and utilize first-order finite-difference donor-cell/upwinding based 
schemes, implemented on structured staggered meshes. 

The schematic presentation of the model hierarchy and relation between semi-implicit algorithms 
is show in the Figure 2. One of the current state of art solution of the problem of the real-time 
control of the two-phase flow in industry is based on the semi-implicit and nearly-implicit algorithms. 
Importantly, the corresponding models were verified and validated for multiple flow regimes and 
incorporate a few decades of extensive numerical experience. 

In this work following the best industrial practice we focus on the development and 
implementation of the semi-implicit and nearly-implicit algorithms for the solution of the fully 
separated non-equilibrium two-phase model of the cryogenic flow. 

Non-hyperbolicity of the Wallis model 

Eigenvalues of the Wallis model 

Before we proceed to the analysis of the solution we note that the set of equations of the Wallis 
model (1) is non-hyperbolic. Two eigenvalues of the system (1) are complex. To see this and to 
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understand the idea behind various techniques introduced in the literature to correct this problem 
we first write system (1 ) in the form 


dF_ dU_ 
dt dx 


( 6 ) 


where 




A 


r g 



ap g u g 


0 


-r g 

(1 -a) Pi 


(1 -a)p,u, 

2 


0 

o 


a P s z ,*-^ s K s -*Ji+T s u ig 

a P s u s 

;U = 

a p, u . 

;C = 

o 

• c o = 

(l-a)pz —t l — r ./ +r u. 

\ ) r g , x wg wg gi i g ig 

(l -a)p,u, 


(l-a)p,uf 


/ \ 


l 1 

a P s e s 




- pcc t -p[au g J x 


Q g w J j + ‘igi-2 +T gK 

(1 -a)p,e,_ 


_(\-a)p,e l u l _ 


pa,-p{(\-a )m ? ) x 


*'* lj A + * u h~ T * ha 


The next step is to linearize this system of equations and write it in the form 


dF_ 

8Q 


C, 


dQ 

dt 


+ 


dU 

dQ 


C 


dQ 

dx 


= c 0 , 


where Q is the vector of primitive variables 

Q-yx p u 


g u i 


e g e - 


O'- 


The eigenvalues of the system (6) are found by solving generalized eigenvalue problem 

/ (A) = det(AA -B) = 0, 


( 8 ) 

( 9 ) 

( 10 ) 


where A = dF/dQ - C t and B = dU/dQ - C x . A Maple code that calculates required Jacobians and 
eigenvalues can be found in Appendix D. The Wallis model has the following set of distinct 
eigenvalues 


X = \ 


u,v, 


ap,+(l-a)p g 


±z^a(l-a) 




ap,+(l-a)p 


( 11 ) 


It is clear that a pair of eigenvalues is complex at any parameters of the model. According to 
Staedke [8]: This indicates that the assumption of exclusively algebraic terms for the interfacial drag 
results in an incomplete formulation for the interfacial momentum coupling. A variety of corrective 
terms is suggested in the literature to mitigate this problem. Below we consider one of such terms, 
which used in the present model. 


Virtual mass term 


There is a common agreement that corrective terms that contain derivatives of the velocity 
or/and pressure are required to complete the formulation of the interfacial momentum transfer. 
However, at present there is no way to deduce these terms completely from the first principles [8], 
A common correction (which is also used in RELAP5 code) is to add so-called the “virtual mass” 
term in the phasic momentum equations accounting for the effect of the local mass displacement in 
the case of a relative acceleration between the two phases. This term has the form (cf Ref [22]) 


M 


V 


Ca(l-a) p m 


d(u g -M z ] 
dt 


dii 

+ u, — - — u 
1 dx 


g 


d u i 

dx 


( 12 ) 
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With an addition of this term the non-hyperbolicity of the Wallis system can be alleviated and 
the eigenvalues of the model take the form 

l[a.p [ u l + (\-a)p g u g ^~ p^i^l-a^v-u^ \fptC 2 («-(l -a)v) 2 -4a(l-a)/o(«-v)(p l)1 CMv + r(M-v)) (13) 

2 [ap, + (1 - a)p g - p m C (1 - a)) 2{ap, + (1 - a)p g - p m C (1 - a)) 




We note that non-hyperbolicity does not necessarily results in the instability, see e.g. [10] for a 
counter-example. However, numerical experience suggests that the introduction of the corrective 
terms is essential for obtaining stable solution. 

A set of equations as it was coded 

In this section, we briefly outline the details of the semi-implicit algorithm as it was coded for the 
chilldown applications. Further details can be found e.g. in [2], [4], [10]. A detailed outline of a 
similar nearly-implicit algorithm is provided in “Part IV: Code Structure” of this report. 

As a first step in deriving this algorithm, we obtain a numerically convenient set of equations for the 
two-phase flow starting from the full set in the form 

(Aa gPg ) t + {AagPgUg}^ = AY g , 


(Aa a p a u a ) „ + ( Aa a p a u l ) .. + Aa a p x = -Aa a p a z mX - T aw l wa - t qi I l + AY a u io , 

(■ Aa s p s E g ), +(^v>A“*), x = A(X p< ~(p Aau s ), x +< iJ™ g + % l i +AY A’ ( 14 ) 

(App a ) „ + {AppiU^ = - AY a , 

(Afipiu i) if . + U/fftu ?)_ + Ap p,* = -App^ - r lw l wl - t u l t - AY a u [b 
(ApEipi) j- + ( ApEip^Ui Appj; {pApiii)^ + "F Quh 

Note that these equations are written for total the energies^ z = e gl + u g,i I 2 E a , = e aml + We 

now consider each pair of equations in turn 


Continuity equations 


The continuity equations have the form 

A (ap g ) t +( A ap g u g ) x =AT g , 

A {PPi ),t + ( A PPi u i ), x = ~ AY g- 

It is convenient to write (15) in the form of the sum and difference of the mass conservation 
equations 

A(ap g +pp,) t + (A(ap g u g + /3p,u,)^ =0, 

A { a Pg “ PPi )j + ( A ( a Ps u g ~ PPi u i )) x = 2AT S • 


(15) 


( 16 ) 


In the derivation of these equations the following jump conditions were used 

r> 

corresponding to the conservation of total mass. 


r =-r 

1 1 1 g> 
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Momenta equations 


The next step is to write sum and difference of the momenta equations in a non-conservative form. 
Expanding momenta equations and subtracting density equations multiplied by corresponding 
velocities we have (in this section ft is used inter-changeably with (1-a)) 


Aa P g U g ,x + \ Aa P g ( U l), 


+ Aap.=-Aap s z, 


' fl 

Jgw 


0 U \U 

g \ g 


- c.u , , \uA + AT 


■(«*-«*) 


+ M„- 


(17) 


A PPi u u + 2 A PP, i u i ) , 


+ A PP,x=~ A PPl Z p 


PL 1 Pi u i p I 


" R R | A ^~ ' g^fli ^/ ) Tf [/ 


where we have introduced velocity difference u R = u a - u t u R =u g -u, term and virtual mass term 
M v = Cafip^iiQ - u; ) M v = aPp m {u g - u , ) t . The latter is a simplified version of the full expression 
(12). 

Next, we notice that it is convenient to cancel cross-sectional area A A whenever possible and 
introduce (see RELAP5-3D [2] for further details) the following definitions of the frictional losses 
terms 




(fl > 

Pg u s 

u g 

l 4 J 

2 


A( 1 a ) Pi u iF w i — I 


( f,L, ) 

Pl U l 

U l 

l 4 J 

2 


Aa Pg F ig ( U g ~ U l) = C gi U R hi, A (! - a )Pl U l F il («; ~ U g) = ~ C li U R \ U R I • 


(18) 


Using these definitions and adding and subtracting two momenta equations we obtain 


a/-’,:", , + + -p("i )_, + ) , + /> . = -P„r,. - a P, U , F ,, -fiP;‘<; F . ~ T , )■ 

Hi Hi 


U -Uj f +- 

g,t i,t 2 


1 — 

\.p s Pi j 


n = u,F , -u F + 

■L ,x l wl g wg 


(19) 


T} U ‘ [U ‘ -P +P - F P’ apjp, 


| + — ^ Mr, 


The following jump conditions were used in derivation of Eqs. (19) 

-Aap g F ig (u g -u,) + AT g (u gi -u g ) + M v - A[Pp l u l F il ( w, - u g ) - AT g (u u -u,)-M v = 0 


corresponding to the fact that all internal interactions must cancel each other according to the 1 st 
Newton’s law. It is further assumed without proof that interfacial drags related to the surface friction, 
mass transfer, and virtual mass cancel each other independently resulting in the following simplified 
set of the jump conditions 

a Pg u g F ig = PPi u i F n ; u gi = u u = u i ■ 

Energy equations 


To simplify and speed up calculations it is convenient [2] to rewrite equations of energy in the 
form that contains only internal energy of liquid and gas. To derive the equations for internal energy 
we first derive equations for the kinetic energy and subtract the kinetic energy equations from the 
equations for the total energy. The equations for kinetic energy are obtained by multiplying 
conservative and non-conservative forms of the momenta equations by 1/2 u and adding them 
together. For example for the kinetic energy of the gas/vapor we obtain 
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11^ ' 

Aap g f 

z j 


+ 


2 A 


u„ 




V 


= y + (~ aA P, f ~ Kj wg ~ T , l < + AT s u i ) 1 


Subtracting this equation from the equation for the total energy and neglecting the interfacial 
transfer of kinetic energy and frictional heating we have 

( Aa P s e s ),, + ( Aa P s e g u g ). A = ~ A P a , ~ p ( Aau s ), A + iJw + ijt + AT A ’ 

( A Ppfi, ) t + ( Afip^u, ) x = -App t - p ( Ap Ul ) x + q lw l wl + q u l t - AT g h a . 

The jump conditions for the energy conservation are 

qji + AT g h ig + q u l t - AT g h u = 0 . 


( 20 ) 


Expanded form of the density and energy equations 

In addition to the set of equations (16), (19), and (20) the first step of both semi-implicit and nearly- 
implicit algorithms requires solution of the expanded form of the coupled sum and difference 
density and energy equations. This set of equations is obtained by expanding time derivatives of 
equations (16) and (20) as follows 

ap s , + Pp u + «, [P g ~ Pi ) ■ + -j( A ( a PP', + PPAi )) r = 0. 


ap g ,t ~ PPi.t +«, [p g + Pi) + -j{ A { a P s u g -f J Pi ll i)) r = 2r V 

(p g e g +p) a ,, + a Pg e g, t + ae gPg,< +\[{ Aa P g e s u s ) >A + P ( Aau g ), A ] = ^„y + % l -j + PgK 

~(pi e i +p)a, + PPi e u +Pe,pL, + 7 [( A PP, e i u i ), x + p{ A P u i ), x If % + A - v A ■ 


( 21 ) 


Equation for the wall temperature 

In the present research the following simplifications are adopted for the analysis of the wall 
temperature. We neglect axial conduction in the wall as compared to the heat transfer through the 
wall on the ambient side and internal flow sides. We also assume that the thickness of the pipe wall 
d w is small enough to estimate pipe wall volume as V w = nDd w , where D D is the pipe diameter. 
The resulting equation for the wall temperature becomes 

Pw c w A w = h gw (T -T w ) + h lw (T l -T w ) + h amb ( T amb - T w ) (22) 


A preview of the algorithm 

The efficiency of the solution in semi-implicit and nearly-implicit algorithms of RELAP5 is archived 
by splitting the operators for conduction in the wall (22) and two-phase flow (16), (19), (20) into 
practically independent problems (see also [4], [23]). In turn, the solution of the latter problem is 
constructed in a two-stage “predictor-corrector” fashion [4]. 

In the first step the expanded mass and energy conservation equations (21) are solved together 
with momenta equations (19). At the second step of the algorithm the unexpanded conservation 
equations (15) and (20) are solved to enforce conservation of the mass and energy. 


10 


At each step of the algorithm the discretized equations are carefully linearized with respect values 
at the new time step to break CFL limitations related to the pressure waves (semi-implicit) and 
material waves (nearby-implicit). 


An outline of the semi-implicit algorithm 


In this section we describe the solution of the Wallis model using semi-implicit algorithm. First, we 
introduce a set of control volumes, and staggered grid. Next, we describe approximation of the 
scalar variables at the interfaces of control volumes and of the velocities at the centers of the 
control volumes required for the integration. Finally, we describe two steps of the solution. The key 
goal of this method is to lift limitations on the CFL number related to the speed of sound. Note, 
however, that the material CFL limitation 


At < 


Ax 


max 



(23) 


still holds in the semi implicit scheme. 

The advantage of the approach proposed in [2] and adopted for the current chilldown model is its 
efficiency and linear scaling of the algorithm with the number of control volumes. These advantages 
are particularly important for the applications to the real time control of cryogenic systems. 

Grid 

To discretize model equations we introduce N control volumes (CV). The scalar values including 
pressure, temperatures, densities, and energies are stored at the centers of the CVs indexed 
L = l,...,N as shown in the Fig. 3 

j- 1 j j+1 



L-2 L-l L L+l 


Figure 3 Control volumes ( L = 1 , . . . , /V) and staggered grid {^j = 1, . . . , N + 1) of the model. 


Upwind values 

The values of the velocities and mass, energy, and would fluxes are calculated at the interfaces 
centered at staggered grid with indexes j = \,...,N + \. The values of scalar variables are 
calculated using upwind scheme as follows 


li 


Pj =}(%-! +< Pl) + \{s uJ +z uJ •5 Pj ,)(%_ 1 -%) = (l + (5. J +\j - s pJ ))^ + (\-{s UJ +z Utj • J PJ ))^- S 


Yj =\Wl-1 + Yl) + \{°uJ +Z U J ^p^Wl^-Yl) 


+ Z UJ- Z PJ 


2 Pg(l),L-\ + P«(Z)J. 


( 24 ) 


1 + ( 


v ; 1 




j* z pj + 


2 A 




P s m,L-i + P S (i),L J 


Yl- 


!-k 




jH* 


j’ z pj + 


2 A 


8(0.H 


Pg(l),L-\ + Pg(l),L 


Yl_ 
2 ' 


Here cp corresponds to the densities while ^corresponds to the void fractions a 9 (i) or energies 
e 3 (i). The values s and zz are defined as follows 

s u j = sign(Uj ); S . = sign(Pj); Z u . = jl if Uj=0, 0 otherwise | ; Z p . = jl if p j= 0, 0 oz/zenW.vej . 

An additional smoothing is introduced for small velocities for scalar values centered at the junctions 

1 when j m > y), m 

< 5 2 (3-2^) w/zen - 7 lim < j m < j lim 
0 vv/;en < - / Hm 




The following notations are used in this smoother £ = (y m + j Um )/( 2 j Um ) , j m = oc jU . + ot, Uj u ltj , and 
= 0.46 m / .V . 


Face centered interpolated values 

In the discretization of the model equations the upwind values introduced above are used in the flux 
terms. For the scalar values entering time derivatives and source terms in momenta equations are 
interpolated between values in neighboring control volumes as follows 

V I ;(i)J = g(i),L-\h-\ + ^g(/yA)/(A- 1 +0- (25) 


Volume-centered velocities 


The volume-centered velocities are required to calculate the kinetic energy in the momenta 
equations. For a simple pipe flow volume-centered velocities can be found using mass 
conservation law 

A a gU),LPg(l),L U gV),L = ^\_^y g (l)jPgU)J U g(l)J + A j+l®g(l)J+lPgU),j+l U g(l)J+l] 


resulting in the following expression 


U g(l),L ~ 


_^jyppjPguij_ 

2 A a n sU)J 

jLJ1 L U g n),LPgO),L 


+ 


^•+i a g(n.y+i^ > g(/),7+i 

^y a g{l\Lpg(l\L 


u g(.l),j + 1 


(26) 


Artificial velocity diffusion 


In the semi-implicit version of the algorithm [2] an additional stabilizing factor is introduced in the 
form of the artificial velocity diffusion 


2 


a gPg Ax 



) 
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The integration of this term on a staggered grid results in the hollowing 
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(apAx) n gj 


8 (A 

u 

u ) 

” / i 

-8 1 Am 

U | 

n 

x \ 

g 

8) 

X 

S) 

L - 1 


expression 


which can be further simplified by writing a finite difference approximation 
1 


F n 

r VG(L),j 


2 V )e 


gU)J 




{ AU g(l )) j+ 1 (^“g(/))| y .j \ U g(l\L-\[[ AU g(l))j (^ W g(/))| y _ i | 


( 27 ) 


First step of the semi-implicit algorithm 

In this section we describe the first step of the algorithm as it is coded. It involves the following sub- 
steps 

• Solve momenta equations for new time velocities in terms of new time pressure 

• Solve expanded equations for the new pressure 

• Find new velocities using new pressure 

• Solve expanded equations for the provisional values of void fractions and energies 

• Find provisional values of the temperatures and interphase mass fluxes 

The main difference with the original version [2] is the use of the incremental values for the new 
scalar variables and velocities, which eliminates the necessity to use special techniques for the 
solution of resulting matrix equations. The notions new and old values refer to the time steps n + 1 
and n correspondingly. 


Discretizing and solving momenta equations 

Integrating the momenta equations (19) over volume of staggered grid 

_R,, d <J 1 + du?' ] Ax, + [jR . ({u 2 g )[-{u 2 g )? - F v ? ; ) + . ((«? )[ - {uf )? - F r ? )) At - 

= - [dpT - dpi* ). At - ( pi - p n L _ x ) At + 


(28) 


-Pljgtej ~ (apj F: gJ (du? + u\j ) - (apj F'j ( du ?' + u ? ) - (du? + u? - duff - u ? ) 


At Ax. 


I \ 2 Y 


t^\ n 

ap 


a p Ji 


p g p, 

((«?);■ 




ap 


— — \ 
Pi~ Px 


PiP> 


(pT'-p? x ^-{f*AKY<j)-f«A^ +u I)- 


(29) 


v y 


p:,jd< -M. du?; -{apt du ? p:Xj -M u? -(apj 


-r" 


H" R ’ 

v ' g,J v 'hj 


R.R. 

v 7 gj v 'hj 

{F,p)“. (du?; + u? - du?; - u ? )} AtAXj 

We note that of the new time step these equations involve only dufj du”! 1 , du ?. du™* 1 , and 
Pl(l~ i)P™a-iv Therefore, at j-th interface the two equations can be presented in the form 

n „«+ 1 . in „«+ 1 . n 


A n 


d <j 

du?' 


=<pr+Kp n L - l +< 
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This equation is solved to express new velocities in terms of new pressure 

7 , 77 + 1 ( An \ ^ n+\ . un n+l . n \ 7 n+1 ( ati\ X n+l . Ln n+ 1 . n \ 

du g J ={ A u \ \ a u p L + b uPL-l+ C u)’ du U ={ A u) 2 [ a uPL + b uPL-l+ C u) • 


( 30 ) 


Note, that in the equations above matrix A" and vectors a n u , b n u , and c” are functions of variables 
at the previous time step. The solution (30) can be substituted into the set of expanded equations to 
find new pressures. 


Solving expanded equations for new pressure 


(31) 


Integrating expanded equations (21 ) over the length of the control volume we obtain for the density 

a lAp n ;l + fadfiC + dafl (pl L - pi ) +y({apA)' t , +] - (^)" , u?J ) + 

y(M L, 'ML <‘) = °* 

a g,APg,L - PijdPu + da g L (p gl + p,j ) + u gj+ 1 - (ap.d)^ J - 

-y(MlM-Ml<Mr:, 


and for the energy 

{plL e lj.+Pl)da n fl+(ap) n gL del}+(ae) n gL dp" g l+— [Aa[pe+ p } ) +[ *C/+t +(da(pe + p'j'j uf' 


M {Kl -T s n ?)s ws+ hi L (f;i -f£)s ig+ ti L hi L ]~, 

-{.Pi.L e ",L + Pl)d a g~L + ( a P) LL del' + ( ae )i, L dp"*L + y (f a (P e+ P)) l/tl w "y+' + (^ a (p e+ 


(32) 


,.n+ 1 

IJ 


\k,l (Kl - Kl ) s* + (zjr 1 - c 1 ) 5, - f ix,l ]y ■ 


V 

To solve these equations we write velocities in the form u H g + 0Xj = u n g(l)j + diy} hj and express 

d u n g + (j)j ' n terms of new pressure using equation (30). Next, we express densities and temperatures 
at the new time step in terms of energies and pressure using Taylor series expansions as follows 


7 n+ 1 _«+l 

dPg,L - Pg,L ~ Pg,l 


dTi' = r: 1 -r, 

g,L g,L g,L 


\ d Pjg,L 
r 8T v 




rap 


fife” 


l de ) "~ g - L ’ 

V y g>i 


dpiL — Pul Pi,l 


f ^ \ n 

dp_ 

V d Pj,,L 


dpT + 


(dp 


de n+x • 




r dT Y 


4T + ^ J <2 ; dT,:i = t ;:; 1 - t;i 


r dT \" 


\ d Pjl,L 


\SeJij. (33) 

f^T de n+1 

[deL 


dpi 1 +1^1 de^. 


The derivatives in these expansions are taken from the NIST tables as will be discussed in more 
details in the Part III of this report that discusses correlations. After these substitutions the system 
of equations can be written as a matrix equation 

1 nin + 1 . un .n + 1 . „nin+\ . jnin+l 


A"dx 


x g,L 


: <^ +1 + + d n x du;y + e ; 


(34) 


for the following unknowns 


r/jc” = , de n yl , de^ , dp n L +1 } 


Below we provide the explicit form of matrixes and vectors in (34) under additional simplifying 
assumption that all the temperatures in the source terms are taken at the old time. The rationale 
behind this assumption is twofold. Firstly, in practical applications the time step for semi-implicit 
algorithm is of the order of 5 ms, which is much longer than the characteristic time of the 
temperature relaxation, since the temperature is the slowest variable of the model. This 
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approximation can be violated for very small void fractions for one of the phases. However, in the 
limit of small void fraction, the whole algorithm breaks down and limiters and smoothers are 
employed to stabilize calculations as will be discussed in the Section related to the phase 
appearance/disappearance. Secondly, the 2 nd step of the algorithm is used to enforce mass and 
energy conservation and the temperature is corrected during this 2 nd step. Ultimately, the final 
approximation will be established after validation of the algorithm and extensive numerical 
experimentation. 

Under these assumptions, the matrix A x has the following form 



-(>■ 


_ n . n 

P g,L + Pl,L 


?d e p) n gJ +Pl L ) 


0 

M\spI 

i ead pp) n g , L 

0 


)( P : L+ (ed eP y lL ) 

-Wu-pl 

((1 -a)ed pP ) n iL 

a lA d 'p) n sJ . 

(>- 

a ^)( d 'p)l 

^ n 

P g,L ~ Pl,L 

a iA d P pf g ^ + ( l - a U)( d P p)[ 


The vectors a" and b " have the form 


■a 


gJ + 1 


Pgj+i 

a .. = _ ({epy^+Pl) 

0 

Ki* 

and the vectors c" and d" have the form 


c, = 


AtA J+ 1 
V, ’ 


b : = 


p g j 

(?p)" gJ + pi) 

0 

Kj 


AtA.. 


a 


g,J 


V, 


P? i 



~PiJ 

0 

At ■ A . , . 


0 

-(i b P)?+Pl) 

• a n 1 ■ 

a, ’P l V ’ 

V L 

d> 


-Pij+i 



L pi.j \ 


■a, 


A t-Aj 


ij 


V T 


Finally, the free vector in the equation (34) is written as follows 

At 


e = 


2T n V 

Z1 g,L V L 


r B h n v -O 

L g,L n wg r L *twg 

_r" k n v -o 

1 g,L n wr L ZZwl 

0 


V, 


+ a x U g,j + 1 + K U g,j + C x U IJ+ 1 + d X U IJ 


To solve equation (34) we notice that on substituting (30) the last row of this equation depends only 


on the new pressure increments dp n L _ x , dp'[ dp[- +l , and dp n L+x dpt H +i ■ Multiplying the last row of eq. 
(34) by inverted matrix A ” we obtain the matrix equation for the pressure in the form 

dp? = a n p dp? + b n p dp? + c n p dp? + d n p . (35) 

This tridiagonal equation can be readily solved using standard algorithm (see Appendix C). 


New velocities and provisional values of densities, void fractions, and energies 

Once new pressure is known the new velocities are obtain by substituting result of solution of (35) 
into (30). Next, using new velocities and pressures we obtain provisional values for densities, void 
fractions, and energies by solving the remaining three equations of (34). Next, we update the 
upwind values of scalar variables using new velocities and pressures. At this stage we also update 
the heat transfer coefficients and interphase mass transfer rate using provisional values for the 
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temperatures and new velocities. This completes the first step of the algorithm that solves momenta 
equations and expanded (non-conservative) form of equations for densities and energies. Since, 
the densities and energies at this step were obtained using solution of the non-conservative form of 
the equations to corrects the errors a second step is required. 

In the following section we describe the calculation of corrected the mass and energies by solving 
unexpanded conservative form of the density and energy equations (15) and (20). 

Second step of the algorithm 


To update densities and energies we integrate equations (15) over the control volume and obtain 
the following equations 


m :; 

Mu 



(36) 


Integrating the energy equations (20) we obtain 


i a P e T g j=( a pei 


g,L 





-a„ 



p n L d&i L + { k S ' L (f: L -f; L )s: g ' L+ v L fjy 


(37) 


(ape)"* = (ape)[ L - a? J+1 f (ep)”^ + plju^A^ ~ f (ep)[. + Pi K/A 


At 

k' 


Pld^" L + (Kl,L (K,L - f",L ) Swl,L ' 


vr n h 

V L L g,L n l,L 


£> 

’V , ' 


In the semi-implicit algorithm these equations are solved explicitly for every control volume. Note 
that the overall implicitness of the solution can be slightly elevated at no cost by taking pressure 
values at the new time step. Note also that various modifications can be justified only after 
extensive numerical experimentation. 

The outlined algorithm is perhaps one of the most efficient approaches to the solution of the two- 
phase Wallis model. It involves inversion of N 4x4 matrices, inversion of N 2x2 matrices, 
solution of tridiagonal matrix equation with NxN matrix, and N x m N x m explicit computations. 
The algorithm scales linearly with the number of control volumes. 

In practice, however, we were not able to exceed time step limit 5 ms. This poses some limitations 
for the real time applications to the autonomous control of cryogenic loading. That is why nearly- 
implicit approach that allows to increase time step several times is considered in the following 
sections. But before we do, let us briefly summarize the steps involved in the semi-implicit 
algorithm. 


Summary of the semi-implicit algorithm 


FIRST STEP (predictions) 

1. Solve sum and difference momenta equations in non-conservative form with respect to pressure corrections 

2. Solve last in the set of expanded equations in with respect to new velocities corrections in terms of pressure. 

3. Substitute new velocities from 2 into solution 1 and solve resulting tridiagonal equation for the pressure 

4. Find new velocities using new pressure obtained at the previous step 

5. Find provisional values for energies, void fractions, and temperatures 

6. Find provisional values of mass exchange rate and heat transfer coefficients using provisional values of 
temperatures obtained in step 5 
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SECOND STEP (corrections): Using conservative forms of the mass and energy conservation 
equations find new values of the 

• densities, 

• void fractions, and 

• energies 


Nearly-implicit algorithm 

The key idea behind the nearly-implicit algorithm is to elevate the implicitness of the scheme to the 
level required to break material CFL number limitations (23). In this respect two main corrections 
were introduced. The kinetic energy in the momenta equations and the conservative forms of the 
density and energy equations are solved implicitly. These corrections require some significant 
changes in the solution strategy (see [2]), which we describe below. 

First step of the nearly-implicit algorithm 

In the first step of the nearly-implicit algorithm we solve the coupled pressure-density-energy 
equations. Although this step is similar to the first step of the semi-implicit algorithm, the non-local 
coupling between pressure and velocity in the momenta equations requires a different order of the 
solution and a different integration scheme for the solution of velocity equations. 


Approximation of the momenta equations for nearly-implicit scheme 


The most significant changes are introduced into the discretized momenta equations (28) and (29) 
and as a result into the solution of the velocity-density-energy coupling equations. The artificial 
diffusion terms (27) are no longer required for stabilization. But the kinetic energy terms are treated 
implicitly. However, in the spirit of maximum computational efficiency these terms are linearized 
with respect to new velocities 


KT=W Ki-Kf * K. 2“«V%'.T +fe )’ = ldu >d +(»R < 38 > 


With these modifications the momenta equations read 
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(40) 
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Note that coupling between velocities and pressure has become nonlocal. Before writing these 
equations in matrix form we express volume centered increments for the velocities du n g + ( ' l)L are 

approximated using eqs. (26). Accordingly, the new velocities in equations (39), (40) are located at 
nodes j - 1, j and j + 1 . Because of the nonlocal character of this coupling the solution of eqs. (39), 
(40) with respect to pressure in terms of velocities becomes impossible. Instead these equations 
are solved as follows. 


Solution of the coupled momenta and expanded density and energy equations 


To solve coupled pressure-density-energy equation in the nearly-implicit scheme we notice that the 
expanded set of equations (34) does not change. Accordingly, the last row of the matrix equation 


a: 


d <l 
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de \ ? 

dpt 


= aidu n g % + b n x du;:; + c n x du^/ +l + d;du \ ;;; + 


n+\ 


n + 1 


n + 1 


can be used to express new pressure in terms of new velocities 

dpi" = (a; )■' (a;du'"„ + bldu’" + cldul", + dldul" A el). 


(41) 


This solution can now be substituted into the momenta equations 
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resulting in the matrix equation for the velocity with A/+1 blocks of size (2x6) of the form (the 
following example is given for N = 3 control volumes) 
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(42) 


The explicit form of the coefficients of this matrix equation is the following. For the 1 st row in (2x6) 
block we have 6 coefficients in the form 
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For the 2 nd row in this block 6 coefficients corresponding to the difference of momenta equations 
are 
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The components of the vector at the right hand side n 2k+l and n 2k have form 
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The new velocities obtained by solving (42) are substituted into (41) to find new pressures. The first 
step of the algorithm is completed by solving the rest of the expanded system of equations for 
density and energy and obtaining provisional values of void fractions, densities, energies, and 
temperatures. This part of the algorithm is equivalent to the semi-implicit scheme. The second step 
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of the algorithm is more involved because of the elevated implicitness of the unexpanded density 
and energy equations as we describe below in more details. 

Second step of the nearly-implicit algorithm 


The second step consists of solution of the un-expanded equations for the density (15) and energy 
(20) conservation using elevated implicitness. Let us consider first gas density equation integrated 
over the control volume. 


d M2 l - ( d M2l + 1 + M2,A u 


•*' A ; * |A< < 43 > 


gJ + 1 


V, 


V, 


Similarly, for the liquid density we have 

d MTl = ( d M2U + M n , J+ 2] u i^ -j2^~( d MTj + M2 )<2—= - r " - At • 


V, 




(44) 


We note that interface values (ap)"^ . are obtained using 2 nd equation in (24). With such a 
substitution the solution of the density equations is reduced to the solution of the tridiagonal matrix. 
The latter is solved using algorithm described in Appendix C. 

For the gas energy equations we obtain 


d M)2l = ( d M)Z + 1 + ML + i + M" gJ+ x — 


At 


V, 


[ d M)2l + M)l + Ml 4^ = ~ p M,L + [k s , L (k L - k L ) 




(45) 


A :i ,At 


Similarly, for the liquid energy we have 

d M)"2 = [ d M)22 + MT lJ+ 1 + M[ J+ , pT ) <j- 

’ L 

( d(ape )j;‘ +(ape)l + (ap) n IJ p?')u%^- = p n L da* giL +{k,,L 


n + 1 j+\ 

7+1 


(46) 


To find primitive variables from the set of new conservative variables 

/ \n + 1 / \n + 1 / \n+l / \n+l] 

M g ,L ’ ’ M) g ,L ’ 


we notice that gas and liquid energies are readily obtained as 

C - {<«*%. 1 ME • <£ - 1 ■ 


To find void fractions and densities are found using 

a i,L ~\ a P)lL Pl,L ’ — a i,L 9 Pg,L ~\ a P)g,L 


n+l / „n+l 
U g,L * 


Phase appearance/disappearance 

Numerical experimentation with semi-implicit and nearly-implicit algorithms reveals that one of the 
most sever stability issues is related to the phase appearance/disappearance. The second critical 
issue in phase-transition is the lack of positivity [25], that frequently appear in the control volumes 
with small mass fraction of one of the phases. In the latter case, even small inaccuracies easily 
lead to negative mass fractions, especially with large time-steps. 
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Phase appearance and disappearance 

According to [4] appearance and disappearance presents very serious numerical challenge for all 
effective field models and no fully satisfactory solutions exist, but several ad hoc numerical 
treatments used in different codes. 

To get further insight into the physical origin of this problem we follow [25] and notice that in the 
limit of vanishing volume fraction of one of the phases, the two phases almost decouple and the 
minority phase obeys a pressureless gas dynamics system. Indeed, consider the reduced model of 
isothermal two-phase flow with the non-hyperbolicity corrected in the form of [26] 

Sa g a l p m (u g -ui) 2 (a g ) x 

( a g P g ), t + ( a g P g u g ), x = 0 

(‘ a 8 P 8 u g ),, + ( a 8 P 8 u l ), x + a g P,* + Sa s a iP m K - u, f (a g ) x = 0 
( a iP g )' t +( a iPi u i), x =° 

{a,p,u, ) t + (a,p,uf ) ^ + a lPx + Sa g a t p m [u g - u, ) 2 (a, ) x = 0 

By introducing a small parameter for the vanishing phase a a = sa aa and using relation 
Ea oa + a t = 1 we have in the limit s -» 0 

{ a 8 oP g ) it +(<* g oP g « g ) 3X =° 

( a g o p g u s ), ( + [<* 8 0 p £ ),, + a g oP, x + sSa g0 a lPm (u g - u, f (a g0 ) x = 0 
(p g )+(p l u t ) x =0 

(. Pi u i ),, + ( Pi u f ) x + P, x ~ s5a g() a l p m (u g - u, ) 2 (a, ) x = 0 

where correction terms also vanish in this limit. As a result, the liquid dynamics is completely 
decoupled and the pressure is determined by the liquid dynamics alone. In this case the gas 
dynamics becomes pressureless, the Jacobian matrix of gas dynamics does not have a complete 
basis of eigenvectors, and the system of equations becomes non-hyperbolic [27], [28]. 

Lack of positivity 

An additional difficulty related to the use of the semi- and nearly-implicit schemes is the lack of 
positivity. The notion of "positively conservative" schemes was introduced in [29], where it was 
shown that while Godunov scheme is positively conservative, i.e. preserve positivity of the density 
and internal energy, no linearized Riemann solver, including the Roe scheme, is positively 
conservative. 

At present no fully satisfactory solutions exist to these problems, but several ad hoc numerical 
treatments are suggested in different codes. In this work, we follow recommendations by Liou [30] 
and adjust temperature, velocity, and density according to the following expression 

tdj = S{x)(l) d +{\-g{x))(l) c -, g(x) = x 2 (2x-3); x= ° d * min , (47) 

T — T 
^max "'min 

where “d” stands for disappearing phase and “c” for conducting phase. The exact values of the 
minimum and maximum void fraction, for which smoothing (47) is applied should be established 
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wining extensive numerical experimentation. The usual values are of the order of x max = 10 -3 and 
x = ltr 7 

The role of the limiters and smoothers described in this section is critical for the code stabilization. 
However, their application is not sufficient to eliminate instabilities. An external loop with time step 
control is required to improve code performance as described in the following section. 


Time step control 


The time step control is critical for reliable performance of the algorithm [2], Its implementation is 
somewhat different in the semi-implicit and nearly-implicit schemes because of the CFL limitations 
related to the material waves in the former scheme. Here we describe common step controls for 
both schemes. The priority is given to the control of the mass and energy conservation. In the 
present version of the algorithm the time-step control follows very closely the recommendations of 
[2]. Mass conservation is controlled in two ways. Firstly, the total density 


rc+1 
P m,L 


*+1 #i+l 
U g,LPg,L 



-a 


n+l 

g,L 



for each control volume is calculated using solution of the unexpanded form of the mass and 
energy conservation equations at the second step of the algorithm (see e.g. eqs. (43)-(46)). These 
total densities are compared to the densities obtained by Taylor expansion (33) 
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in the first step of the algorithm. The maximum error is found as follows 


Err = max 


f «+l ~ «+l 

Pm,L ~ Pm,L 


P, 


n+l 

m,L 


Secondly, the error for the total mass in the pipe is found as follows 



If either of these errors is larger than 10 -z the time step is halved. 

In addition, the total mass “in” and “out” of the pipe through the input/output and damp valves is 
continuously monitored and compared with the total changes of the mass within the pipe. If the 
difference exceed a few percent of the total mass, the integration is terminated. 

We also monitor the mass addition/removal for each phase in every control volume. If mass 
addition/removal exceeds half of the mass currently present in the control volume the time step is 
halved. 

Despite this control the values of the key thermodynamic variables may occasionally be found to lie 
outside predefined range. If this happens the time step is halved. The reduction of the time step 
continues until At reached the minimum predefined value (usually 10 mks). If all the field values are 
found to be inside the range, and both mass errors are smaller than 10 -3 , the time step is doubled 
until At reached the maximum predefined value (usually 10 ms). 


Conclusions 


Using limiters and smoothers to control of the phase appearance/disappearance in conjunction with 
external loop for the time-step control a reliable and stable performance of the semi-implicit and 
nearly-implicit algorithms is achieved. 
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Note that a significant part of the algorithm related to the suitable choice of the source terms and 
the discussion of the verification of the algorithm performance will be discussed in separate parts of 
this report. The reason is that the functional form of the source terms depend ssubstantially on the 
application of the algorithm. Furthermore, currently the exact form of the source terms is not 
available/known for the application to the real time control of cryogenic loading. Accordingly, in the 
Part III of the report we will discuss a few possible approaches to the solution of this problem. In 
addition, we will provide example of a full set of complete correlations for the flow patterns, friction 
losses and heat transfer that were developed for two industrial applications. 

Appendices 

Appendix A. Equations of the two-phase flow 

Various approaches can be utilized in the derivation of the two-phase flow equations (see e.g. [8]- 
[1 1], [18], [19]). Here we follow mainly derivation proposed by H. Staedtke [8], which in turn follows 
closely derivation [20] based on the averaging method [21]. We note that this derivation does not 
provide first principles proof of the two-phase flow equations. However, it highlight the key 
properties of the flow that have to be understood for such proof to become possible. The method 
[20], [21] can be shortly summarized as follows (repeated from [31]) 

. . fl, if x belongs to phase k, 

• Introduce phase indicator % k \ x,t\ = < Note that the gradient of the 

[ 0, otherwise. 

distribution function is zero everywhere except the interface. Then the unit vector at the interface is defined 
as n'f = — (V y k /\V y k | j . The Vy fc is assumed to have properties of the delta-function. 


• then Gauss theorem reads I//') = (V (% k ys ~ (ig V ' % k 

• and Leibniz rule reads (x k ¥ t ) = {{x k ¥ ) t ) ~ (x k ) t ) 

• taking into account {ad hoc no surface break up or coalition) that for an observer moving with the interface 
Zfc does not change in time ( Xk ), + u o^Xk = 0 • 

• The mass conservation after averaging takes the form 

{XkPj ) + (ZtV (pu)) = (( XkP l ) - (p (Xk X ) + ( v (x k P u )) - {p uV {Xk)) = Q or 

({XkP), t ) + {v{XkP u )) = (p{Xk), l ) + (p uV {Xk)) and takin 9 int0 account equation for we 

have ([XkPk f } + ( V ( XkPk u k )) = (Pk ( u k ~ ll o ) V ( X k )) • Tn derivation we take into account that %h 

is nonzero only when p = . Note that the term on the right hand side is the volumetric source term for the 

inter-phase mass flow into phase k. 

• Similarly for the momentum equation we have 

(XkPk U k ) + V {XkPk U k ) + v {XkP) = {p^Xk ) + { Pk U k i u k - W 0 ) v Xk ) 

• And for the total energy we have 

(x k Pk E k } + v {XkPk u k ( E k + p/pk )) = (p u k V Xk } + ( Pk E k ( u k - u o ) V Xk ) 

Note that this derivation reproduces non-hyperbolic character of the model. Note also terms that 
contain Vj* and correspond to the interfacial contributions. Specifically, terms proportional to the 



To shorten notations introduce volume average as — — 
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p k (u k ~ u a)V(Xk) correspond to the terms in the model (14) that are proportional to the interfacial 
mass flow r a . 


Appendix B. Single-phase flow equations. 

To contrast a formal derivation above to a simplified derivation of the single- and two-phase flow 
in a pipe we follow here a discussion in [19]. For simplicity, we consider rectangular pipe with 
stratified flow shown in the Fig. 4. Let us assume initially a single-phase flow when the colored part 
of the flow represents the entire control volume with dimensions Ax 1 ,Ax 2 , and Ax 3 and flow direction 
along x x . 



Figure 4 Rectangular pipe with stratified flow. 


In this case mass conservation for a given control volume reads 

AXjAx 2 Ax 3 = j\ (xj)Ax 2 Ax 3 - j\ (xj + Axj)Ax 2 Ax 3 = ( pu 


-(pu) 




Ax 2 Ax 3 . 


Let us consider now the case of two-phase flow with interfacial mass flux. Because of the mass 
exchange, the volume of liquid will also change as shown by the red dashed lines. In this case Ax 3 
becomes function of time and the mass conservation for liquid in a given control volume takes the 
form 


2-(/?(r + Ar)(Ax 3 +u 0 Ar)-/?(r)Ax 3 )AXjAx 2 = /, (xj)Ax 2 — j \ (xj + Axj)Ax 2 Ax 3 +r ; AF 

where u 0 Atis the change of the Ax 3 , u 0 is the velocity of the interface surface, and the interfacial 
mass flow is given by 

T ; AF = pu 0 AxjAx 2 . 


Here Ax t Ax 2 is the interface area and r ; is normalized to the control volume. 

Appendix C. Standard tridiagonal solver 

void solve_tridiagonal_in_place_destructive (double x[], const size_t N, const 
double a[], const double b[] , double c[]) { 

/* unsigned integer of same size as pointer */ 
size_t in; 

/* 

solves Ax = v where A is a tridiagonal matrix consisting of vectors a, b, c 
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note that contents of input vector c will be modified, making this a one-time-use function 
x[] - initially contains the input vector v, and returns the solution x. indexed from [0, ..., N - 1] 

N - number of equations 

a[] - subdiagonal (means it is the diagonal below the main diagonal) — indexed from [1, N - 1] 

b[] - the main diagonal, indexed from [0, ..., N - 1] 

c[] - superdiagonal (means it is the diagonal above the main diagonal) — indexed from [0, ..., N - 2] 
*/ 

c [ 0 ] = c [ 0 ] / b [ 0 ] ; 
x [ 0 ] = x [ 0 ] / b [ 0 ] ; 

/* loop from 1 to N - 1 inclusive */ 
for (in = 1; in < N; in++) { 

double m = 1.0 / (b[in] - a[in] * c[in - 1]); 
c[in] = c[in] * m; 

x[in] = (x[in] - a[in] * x[in - 1]) * m; 

} 

/* loop from N - 2 to 0 inclusive, safely testing loop end condition */ 
for (in = N - 1/ in-- >0/ ) 

x[in] = x[in] - c[in] * x[in + 1] ; 

} 


Appendix D. Hyperbolic analysis of the model equations 

Below we use vector V where we substituted alpha with r only to avoid differentiation of alpha 
b := 'b'; 

U := [alpha, u, h, v, p, H]; 

V := [r, u, h, v, p, H]; 

a := Vector([alpha*rho, alpha*rho*u, alpha*rho*h, (1-alpha)*r, (1-alpha)*r*v, (1-alpha)*r*H]); 
b := Vector([alpha*rho*u, alpha*rho*u A 2, alpha*rho*h*u, (1-alpha)*r*v, (1-alpha)*r*v A 2, (1- 
alpha)*v*r*H]); 

c := Vector([0, 0, -alpha*p, 0, 0, -(1-alpha)*p]); 
d := Vector([0, alpha*p, 0, 0, (1-alpha)*p, 0]); 

A := Jacobian(a, U); 

S := Jacobian(c, V); 

M := A-S; 

B := Jacobian(b, U); 

Y := Jacobian(d, V); 

N := B-Y; 

Ev := simplify(Eigenvalues(M, N)); 
egv := Ev[[3, 4]]; 
egv := Ev[[1,2]]; 
egv := Ev[[5, 6]] 

Cm := 'Cm'; 
x := 'x'; 

Cm := R*x*(1-x)*rm*(u-v); 

a := Vector([alpha*rho, alpha*rho*u+Cm, alpha*rho*h, (1-alpha)*r, (1-alpha)*r*v-Cm, (1- 
alpha)*r*H]); 

b := Vector([alpha*rho*u, alpha*rho*u A 2, alpha*rho*h*u, (1-alpha)*r*v, (1-alpha)*r*v A 2, (1- 
alpha)*v*r*H]); 

c := Vector([0, 0, -alpha*p, 0, 0, -(1-alpha)*p]); 
d := Vector([0, alpha*p, 0, 0, (1-alpha)*p, 0]); 

A := Jacobian(a, U); 
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S := Jacobian(c, V); 

M := A-S; 

B := Jacobian(b, U); 

Y := Jacobian(d, V); 

N := B-Y; 
x := alpha; 

Ev := simplify(Eigenvalues(M, N)); 
egv := Ev[[3, 4]]; 
egv := Ev[[1,2]]; 
egv := Ev[[5, 6]] 

S := 'S'; 

S := R A 2*alpha A 2*rm A 2*u A 2-2*R A 2*alpha A 2*rm A 2*u*v+R A 2*alpha A 2*rm A 2*v A 2- 
2*R A 2*alpha*rm A 2*u A 2+2*R A 2*alpha*rm A 2*u*v- 

4*R*alpha A 2*r*rm*u*v+4*R*alpha A 2*r*rm*v A 2+4*R*alpha A 2*rho*rm*u A 2- 

4*R*alpha A 2*rho*rm*u*v+R A 2*rm A 2*u A 2+4*R*alpha*r*rm*u*v-4*R*alpha*r*rm*v A 2- 

4*R*alpha*rho*rm*u A 2+4*R*alpha*rho*rm*u*v+4*alpha A 2*r*rho*u A 2- 

8*alpha A 2*r*rho*u*v+4*alpha A 2*r*rho*v A 2-4*alpha*r*rho*u A 2+8*alpha*r*rho*u*v-4*alpha*r*rho*v A 2; 
si := collect(R A 2*alpha A 2*rm A 2*u A 2-2*R A 2*alpha A 2*rm A 2*u*v+R A 2*alpha A 2*rm A 2*v A 2, [R A 2, 
alpha A 2, rm A 2], 'recursive'); 

s2 := collect(-4*alpha*r*rho*u A 2+8*alpha*r*rho*u*v-4*alpha*r*rho*v A 2, [r, alpha, rho], 'recursive'); 
s3 := collect(4*alpha A 2*r*rho*u A 2-8*alpha A 2*r*rho*u*v+4*alpha A 2*r*rho*v A 2, [r, alpha A 2, rho], 
'recursive'); 

s31 := ((-alpha A 2+2*alpha-1 )*4)*R*rho*rm*u A 2-4*R*rho*rm*u(alpha-1 ) A 4; 
s5 := -2*R A 2*alpha*rm A 2*v A 2-4*R*alpha A 2*r*rm*v A 2+R A 2*rm A 2*v A 2; 
s6 := 2*R A 2*alpha*rm A 2*u*v+4*R*alpha A 2*r*rm*u*v+4*R*alpha A 2*rho*rm*u*v- 
8*R*alpha*rho*rm*u*v; 
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Abbreviations 

AUSM — Advection Upstream Splitting Method schemes for hyperbolic systems of conservation laws which do not require any 

analytical calculation of the Jacobian 

RELAP — Reactor Excursion and Leak Analysis Program 

TRAC — Transient Reactor Analysis Code 

TRACE — TRAC/RELAP Advanced Computational Engine 
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